%% Treynor-Black portfolio optimization
% Computes Treynor-Black optimal portfolio weights
%% Rights and disclaimer
% Author: Benjamin J. J. Voigt (bvoigt@gmail.com)
% Published: Oct 2012
% Version: 0.0.2
% License: BSD 2-Clause, free to use and copy retaining original author reference
% 
% Disclaimer: USE AT YOUR OWN RISK; NO INVESTMENT RECOMMENDATION IMPLIED;
% NOT ASSOCIATED WITH PAST, CURRENT OR FUTURE EMPLOYER OR ASSOCIATES OF THE
% AUTHOR; DO NOT USE THIS CODE TO INVEST YOUR OWN MONEY OR TO 
% ENCOURAGE OTHERS TO INVEST MONEY;
%% Syntax
% aaa
%% Arguments
% bbb
%% Examples
% ccc
%% Required toolboxes
% Financial Toolbox, Data Feed Toolbox, Statistics
% Toolbox
%% Purpose and notes
% Purpose: May this code evolve and be enhanced by many smarter minds to
% reflect all the beauty of estimations we dare to dream of.
%
% Notes: --
%
% Known limitations:
% 1)The index date array is used to align individual stocks, actual trade
% dates may NOT match. The reason is data problems in some sources
% regarding dates and the low importance of dates for the statistical
% calculations performed by the function
% 2)Annualization from monthly to 1 year yields estimates for a 1 year
% investment horizon and holding period only
% 3)To improve readability many for-loops are used which limits application
% to a larger number of portfolio holdings without optimization
% 4)The sequence of symbols cannot change during the functions runtime
% or else the weight calculation may fail or result in wrong results
% 5)No constraints on wights may result in abnormal high long/short
% poistions recommended, i.e, in excess of 100%
%

function [g_symbol_optimal_weights g_active_port_alpha g_active_port_beta...
    g_active_port_pct_of_total g_passive_port_pct_of_total, varargout]...
    = TreynorBlack(price_fints, return_fints, symbol_cell, index_cell, ...
    risk_free_rate, symbol_alpha_forecast_mat,...
    symbol_alpha_forecast_acuracy_cell, frequency_str)

% Global, function wide variables
g_symbol_cell = symbol_cell;
g_symbol_alpha_estimates_mat = symbol_alpha_forecast_mat;
g_rfr = risk_free_rate;
g_index_cell = index_cell;
g_price_fints = price_fints;
g_return_fints = return_fints;
g_reg_results = (zeros(length(g_symbol_cell),19));
g_symbol_riks_para = (zeros(length(g_symbol_cell),4));
g_index_risk_para = (zeros(1,4));
g_frequency_str = frequency_str;

g_symbol_optimal_weights = 0;
g_active_port_alpha = 0;
g_active_port_beta = 0;
g_active_port_pct_of_total = 0;
g_passive_port_pct_of_total =  0;

% Local, dummy, loop or sub-section specific variables
l_data_cell = (zeros(length(g_price_fints),2));
l_reg_mat = (zeros(length(g_return_fints),2));
l_whichstats ={'R', 'rsquare', 'adjrsquare', 'fstat', 'tstat', 'mse'};
l_reg_stats = (0);
l_active_weights_mat = zeros(length(g_symbol_cell),1);
l_active_weights_squared_mat = zeros(length(g_symbol_cell),1);
l_std_alpha_mat = zeros(length(g_symbol_cell),1);
l_active_port_alpha = 0;
l_active_port_resid_var = 0;
l_active_port_beta = 0;
l_w_0 = 0;
l_active_port_pct_of_total = 0;
l_sum_of_std_alpha = 0;
% TODO: provide more suitable varargout definition, this will triger
% varargout if-condition all of the time. Also enriching output with
% descriptive cell headers may help.
varargout = {'1' '2' '3' '4' '5' '6' '7'};

% Run regressions
for i=1:length(g_symbol_cell)
    l_reg_mat = fts2mat([g_return_fints.(strrep(cell2mat(g_index_cell(1)),'^','')) ...
        g_return_fints.(strrep(cell2mat(g_symbol_cell(i)),'-',''))]);
    l_reg_stats = regstats(l_reg_mat(1:end,2), l_reg_mat(1:end,1), 'linear', l_whichstats);
    % 1)rsq,2)adj_rsq,3)SE,4)av_DF_R,5)av_DF_E,6)av_SS_R,7)av_SS_E,8)av_MS_R,...
    % 9)av_MS_E,10)av_F,11)av_F_p,12)alpha,13)beta,14)a_se,15)a_t,16)a_p,17)b_se,18)b_t,19)b_p
    g_reg_results(i,1:19) = [l_reg_stats.rsquare l_reg_stats.adjrsquare...
        (l_reg_stats.fstat.sse/l_reg_stats.fstat.dfe)^.5 l_reg_stats.fstat.dfr...
        l_reg_stats.fstat.dfe l_reg_stats.fstat.ssr l_reg_stats.fstat.sse...
        l_reg_stats.fstat.ssr/l_reg_stats.fstat.dfr ...
        l_reg_stats.fstat.sse/l_reg_stats.fstat.dfe ...
        l_reg_stats.fstat.f l_reg_stats.fstat.pval ...
        l_reg_stats.tstat.beta(1,1) l_reg_stats.tstat.beta(2,1)...
        l_reg_stats.tstat.se(1,1) l_reg_stats.tstat.t(1,1) l_reg_stats.tstat.pval(1,1)...
        l_reg_stats.tstat.se(2,1) l_reg_stats.tstat.t(2,1) l_reg_stats.tstat.pval(2,1)];
        
        
end

% Annulization of regression results and risk parameter collection
% sdt_excess_return = SDT(i)*SQR(12), beta, sdt_sys = beta*sdt_INDEX_excess_return, 
% sdt_resid = SE(i-stock)*SRQ(12)

switch g_frequency_str
    case 'd'
        g_index_risk_para(1,1:4) = [...
            std(fts2mat(g_return_fints.(strrep(cell2mat(g_index_cell(1)),'^',''))))*(220^(1/2)) ...
            1 ...
            1*(std(fts2mat(g_return_fints.(strrep(cell2mat(g_index_cell(1)),'^',''))))*(220^(1/2))) ...
            0
            ]
        for i=1:length(g_symbol_cell)
            g_symbol_riks_para(i,1:4) = [...
                std(fts2mat(g_return_fints.(strrep(cell2mat(g_symbol_cell(i)),'-',''))))*(220^(1/2))...
                g_reg_results(i,13)...
                g_reg_results(i,13)*g_index_risk_para(1,1)...
                g_reg_results(i,3)*(220^(1/2))
                ]
        end
    case 'm'
        g_index_risk_para(1,1:4) = [...
            std(fts2mat(g_return_fints.(strrep(cell2mat(g_index_cell(1)),'^',''))))*(12^(1/2)) ...
            1 ...
            1*(std(fts2mat(g_return_fints.(strrep(cell2mat(g_index_cell(1)),'^',''))))*(12^(1/2))) ...
            0
            ]
        for i=1:length(g_symbol_cell)
            g_symbol_riks_para(i,1:4) = [...
                std(fts2mat(g_return_fints.(strrep(cell2mat(g_symbol_cell(i)),'-',''))))*(12^(1/2))...
                g_reg_results(i,13)...
                g_reg_results(i,13)*g_index_risk_para(1,1)...
                g_reg_results(i,3)*(12^(1/2))
                ]
        end
    case 'a'
        g_index_risk_para(1,1:4) = [...
            std(fts2mat(g_return_fints.(strrep(cell2mat(g_index_cell(1)),'^','')))) ...
            1 ...
            1*(std(fts2mat(g_return_fints.(strrep(cell2mat(g_index_cell(1)),'^',''))))) ...
            0
            ]
        for i=1:length(g_symbol_cell)
            g_symbol_riks_para(i,1:4) = [...
                std(fts2mat(g_return_fints.(strrep(cell2mat(g_symbol_cell(i)),'-',''))))...
                g_reg_results(i,13)...
                g_reg_results(i,13)*g_index_risk_para(1,1)...
                g_reg_results(i,3)
                ]
        end
    otherwise
        % TODO: error
        ;
end
        
        


% Calculate weight support values
% std_resid^, alpha/std_resid^2
for i=1:length(g_symbol_cell)
    l_std_alpha_mat(i,1:1) = [g_symbol_alpha_estimates_mat(i)/(g_symbol_riks_para(i,4)^2)];
end
l_sum_of_std_alpha = sum(l_std_alpha_mat);

% Calculate active weights
% w(i), w(i)^2
for i=1:length(g_symbol_cell)
    l_active_weights_mat(i,1:1) = [(l_std_alpha_mat(i,1)/l_sum_of_std_alpha)];
    l_active_weights_squared_mat(i,1:1) = [(l_std_alpha_mat(i,1)/l_sum_of_std_alpha)^2];
end

% Calculate total weights
% TODO: re-check if beta is correct
% TODO: insert variable constraint on wA (= l_active_port_pct_of_total)<=1
%       using factor limit of short in passive, e.g., 0.4691163/(5.7937*(1-factor))
% Note: the exactly composition of input matrices is important because of
% the matrix transpose operator in some of the equations below.
l_active_port_alpha = g_symbol_alpha_estimates_mat*l_active_weights_mat;
l_active_port_resid_var = (g_symbol_riks_para(1:end,4).^2)'*l_active_weights_squared_mat;
l_w_0 = (l_active_port_alpha/l_active_port_resid_var)/(g_rfr/...
    ((std(fts2mat(g_return_fints.(strrep(cell2mat(g_index_cell(1)),'^',''))))*(12^(1/2)))^2));
l_active_port_beta = g_reg_results(1:end,13)'*l_active_weights_mat;
l_active_port_pct_of_total = l_w_0/(1+(1-l_active_port_beta)*l_w_0);

% Assign primary result variables
g_symbol_optimal_weights = l_active_weights_mat(1:end,1:1)*l_active_port_pct_of_total;
g_active_port_alpha = l_active_port_alpha;
g_active_port_beta = l_active_port_beta;
g_active_port_pct_of_total = l_active_port_pct_of_total;
g_passive_port_pct_of_total = 1-l_active_port_pct_of_total;

% Assign secondary, optional, output variables
if(length(varargout)>0)
    switch length(varargout)
        case 1
            varargout(1) = {g_reg_results}
        case 2
            varargout(1) = {g_symbol_riks_para}
            varargout(2) = {g_index_risk_para}
        case 3
            varargout(1) = {g_symbol_riks_para}
            varargout(2) = {g_index_risk_para}
            varargout(3) = {g_reg_results}
        case 7
            varargout(1) = {g_symbol_riks_para}
            varargout(2) = {g_index_risk_para}
            varargout(3) = {g_reg_results}
            varargout(4) = {l_std_alpha_mat}
            varargout(5) = {l_active_weights_mat}
            varargout(6) = {l_active_port_resid_var}
            varargout(7) = {l_w_0}
        otherwise
            % TODO: Error code/message
            ;
    end

end







